# Load necessary libraries
if(!require(pacman)) install.packages("pacman") 
pacman::p_load(
  tidyverse,    # for data manipulation
  data.table,   # for reading data tables
  skimr,        # for checking missingness of data
  gt,           # for creating tables
  rmdformats,   # for changing format of the html document
  gtExtras,     # theming package for gt tables
  lubridate,    # for date manipulation
  sf,           # for creation of spatial visualization
  naniar,       # for visualizing missingness of data
  stringr,      # for string manipulation
  geofacet,     # for data visualization of plots in US shapes
  remotes,      # for installing packages from github
  leaflet,      # for mapping casualty data on map
  scales,       # for scaling plot axes
  patchwork,    # for combining multiple plots
  htmltools     # for HTML manipulation
)

ggplot2::theme_set(ggthemes::theme_few())

## Reading the gun data

gun_data <- fread("Data/gunViolenceGeo.csv ")

glimpse(gun_data)

## Reading the us county data

us_census_data<-fread("https://raw.githubusercontent.com/dilernia/STA418-518/main/Data/census_data_county_2009-2021.csv")

glimpse(us_census_data)

## Checking the missingness of the data using skim from the skimr package
## First converting the date 

gun_data$date<-lubridate::as_date(gun_data$date)

skim(gun_data)

# using naniar to visualize the data missingness

gg_miss_var(gun_data, show_pct = TRUE)

# Using gg_miss_fct to visualize data

gun_data |>
  dplyr::select(state, n_killed, n_injured, participant_type,participant_gender,participant_age_group) |>
  gg_miss_fct(fct = state)

## Dropping the columns

gun_data<-gun_data |>
  dplyr::select(-c(participant_relationship, 
                   location_description,
                   participant_name,
                   state_house_district,
                   state_senate_district
                   ))

duplicated_rows_gun_data<-gun_data[duplicated(gun_data),]
duplicated_rows_gun_data

duplicated_rows_census_data<-us_census_data[duplicated(us_census_data),]
duplicated_rows_census_data

# Converting the year column to date

us_census_data$year<-lubridate::make_date(us_census_data$year)
us_census_data$year<-lubridate::year(us_census_data$year)

gg_miss_var(us_census_data, show_pct = TRUE)

grouping_states<-gun_data|>
  group_by(state,year = year(date))|>
  summarize(total_injured = sum(n_injured),
            total_killed = sum(n_killed))

population_group<-us_census_data |>
  group_by(state,year)|>
  summarize(population = sum(population))

## Adding the population

grouping_states<-grouping_states|>
  left_join(population_group)

## Calculating the rates per 100,000

grouping_states<-grouping_states|>
  mutate(
    injured_per_100k = total_injured/population * 100000,
    killed_per_100k = total_killed/population * 100000)

## Getting the mean per state

mean_states<-grouping_states|>
  group_by(state)|>
  summarize(mean_injured_rate = mean(injured_per_100k),
            mean_kill_rate = mean(killed_per_100k))

## Creating a table of the states with highest casualty rates

mean_states|>
  arrange(desc(mean_injured_rate))|>
  slice_head(n=15)|>
  gt()|>
  cols_label(
    state = "State",
    mean_injured_rate = "Injured per 100k",
    mean_kill_rate = "Killed per 100k"
  )|>
  tab_header(
    title = "Top 15 states with highest injuries/deaths per state between 2013 and 2018"
  )|>
  gtExtras::gt_theme_538()

## Dumbbell plot of states' casualty rates 

mean_states |>
  ggplot(aes(x = mean_injured_rate, 
             y = reorder(state, mean_injured_rate))) +
  geom_segment(aes(x = 0, 
                   xend = mean_injured_rate, 
                   yend = state)) +
  geom_point(aes(x = mean_injured_rate, 
                 color = "People Injured"), 
             size = 4, 
             alpha = 0.5) +
  geom_point(aes(x = mean_kill_rate, 
                 color = "People Killed"), 
             size = 4, 
             alpha = 0.5) +
  scale_alpha_continuous(range = c(0, 1)) +
  scale_color_manual(values = c("People Injured" = "#192bc2", "People Killed" = "#ff0000"),
                     guide = guide_legend(title = "Guide")) +
  labs(
    x = "Mean Injured Rate",
    y = "State",
    caption = "Source: Gun Violence Archive’s website",
    title = "Gun Violence by State per 100k people between 2013 and 2018"
  ) 

## Grouping the casualty rates by month

grouping_month<-gun_data|>
  group_by(month = lubridate::floor_date(date,"month"))|>
  summarize(total_injured = sum(n_injured), total_killed = sum(n_killed))

## Creating an area plot of the casualties

grouping_month |>
  ggplot(aes(x = month)) +  
  geom_area(aes(y = total_injured, fill = "People Injured"), stat = "identity") +
  geom_area(aes(y = total_killed, fill = "People Killed"), stat = "identity") +
  scale_y_continuous(expand = expansion(mult=c(0,0.1)))+
  scale_fill_viridis_d(option = "plasma")+
  labs(
    x = "Month/Year",
    y = "Frequency",
    caption = "Source: Gun Violence Archive’s website",
    title = "Number of gun related injuries/death per month",
    fill = "Guide"
  ) 

## Checking the days in which most casualties occurred


grouping_day<-gun_data|>
  group_by(day = lubridate::wday(date, abbr = FALSE, label = TRUE))|>
  summarize(total_injured = sum(n_injured), total_killed = sum(n_killed)) |>
  arrange(desc(total_injured))|>
  pivot_longer(cols = c(total_injured, total_killed),
               names_to = "Type",
               values_to = "Count")

grouping_day|>
  ggplot(aes(x = day, y = Count, fill = Type)) +  
  geom_bar(stat = "identity", position = "dodge") +
  scale_y_continuous(expand = expansion(mult = c(0, 0.1))) +
  scale_fill_viridis_d(option = "plasma") +
  labs(
    x = "Day",
    y = "Frequency",
    caption = "Source: Gun Violence Archive’s website",
    title = "Number of gun-related injuries/death per day",
    fill = "Guide"
  ) 

## Creating casualty rate visualization for seasons


gun_data <- gun_data |>
  mutate(months = as.integer(format(date, "%m")), 
         season = cut(months, 
                      breaks = c(0, 2, 5, 8, 11, 12), 
                      labels = c("winter", "spring", "summer", "fall", "winter"),
                      include.lowest = TRUE)
  )

grouping_seasons<-gun_data|>
  group_by(season)|>
  summarize("Total injured" = sum(n_injured), "Total killed" = sum(n_killed))|>
  arrange(desc("Total injured"))|>
  pivot_longer(cols = c("Total injured", "Total killed"),
               names_to = "Type",
               values_to = "Count")

grouping_seasons|>
  ggplot(aes(x = season, y = Count, fill = Type)) +  
  geom_bar(stat = "identity", position = "stack") +
  scale_y_continuous(expand = expansion(mult = c(0, 0.1))) +
  scale_fill_viridis_d(option = "plasma") +
  labs(
    x = "Season",
    y = "Frequency",
    caption = "Source: Gun Violence Archive’s website",
    title = "Number of gun-related injuries/death per season",
    fill = "Guide"
  ) 

## Creating a pie chart to display the age group, gender, harmed and unharmed


total_harmed<-sum(gun_data$n_killed, gun_data$n_injured,na.rm = TRUE)
total_unharmed<-sum(gun_data$n_unharmed,na.rm = TRUE)

total_children<-sum(gun_data$participant_child,na.rm = TRUE)
total_teen<-sum(gun_data$participant_teen,na.rm = TRUE)
total_adult<-sum(gun_data$participant_adult,na.rm = TRUE)

total_males<-sum(gun_data$participant_males,na.rm = TRUE)
total_females<-sum(gun_data$participant_females,na.rm = TRUE)

# Creating pie chart for harmed and unharmed
harmed_unharmed<-data.frame(status = c("harmed","unharmed"),
                            count = c(total_harmed, total_unharmed))

chart_harmed<-ggplot(data = harmed_unharmed,
       aes (x = "",
             y = count,
            fill = status))+
  geom_bar(stat="identity", width=1, color="white") +
  coord_polar("y", start=0) +
  
  labs(title = "Pie chart of harmed vs unharmed") +
  scale_y_continuous(labels = comma)+
  scale_fill_viridis_d(option = "plasma")+
  theme_void() 

# Creating pie chart for gender 
gender<-data.frame(status = c("male","female"),
                            count = c(total_males,total_females))

chart_gender<-ggplot(data = gender,
       aes (x = "",
             y = count,
            fill = status))+
  geom_bar(stat="identity", width=1, color="white") +
  coord_polar("y", start=0) +
  
  labs(title = "Pie chart of gender") +
  scale_y_continuous(labels = comma)+
  scale_fill_viridis_d(option = "plasma")+
  theme_void() 

# Creating pie chart for age group
age_group<-data.frame(status = c("adult","teen","child"),
                            count = c(total_adult, total_teen,total_children))

chart_age<-ggplot(data = age_group,
       aes (x = "",
             y = count,
            fill = status))+
  geom_bar(stat="identity", width=1, color="white") +
  coord_polar("y", start=0) +
  
  labs(title = "Pie chart of age group") +
  scale_y_continuous(labels = comma)+
  scale_fill_viridis_d(option = "plasma")+
  theme_void() 

## Combining  the graphs together
chart_age + chart_gender  + chart_harmed 

## Visualization of monthly casualty rates on states plot


time_series_state<-gun_data|>
  group_by(state,date)|>
  summarize(total_casualties = sum(n_injured+n_killed))

time_series_state<-time_series_state |> mutate(year = year(date)) |>
  left_join(population_group, by = c("year","state"))|>
  mutate(casualties_per_100k = (total_casualties/population) * 100000)

sea_palette <- c("#005a70", "#0088a1", "#27a4f8", "#00d0d0", "#ffffff")

time_series_state|>
  ggplot(aes(x = date, y = casualties_per_100k))+
  ylim(0,1)+
  labs(
    title = "Number of gun-related casualties between 2013 and 2018",
  )+
  geom_line(color = sea_palette[3]) + 
  facet_geo(~state)+
  theme(
    plot.title = element_text(color = sea_palette[1], size = 14, face = "bold", margin = margin(b = 10)),
    axis.title.y = element_text(color = sea_palette[2], size = 12, face = "bold"),
    axis.text = element_text(color = sea_palette[4], size = 10),
    panel.background = element_rect(fill = sea_palette[5]),
    strip.background = element_rect(fill = sea_palette[3], color = "white"),
    strip.text = element_text(color = "white", size = 10),
    legend.position = "none" 
  )

### Making a leaflet plot

 if(!file.exists("cb_2018_us_state_500k.zip")) {
  download.file(url = "https://www2.census.gov/geo/tiger/GENZ2018/shp/cb_2018_us_state_500k.zip",
              destfile = "cb_2018_us_state_500k.zip")
}

# # Create directory for geospatial files
if(!dir.exists("GeoFiles")) {
 dir.create("GeoFiles")
 }

 # Unzipping files
 utils::unzip("cb_2018_us_state_500k.zip",
              exdir = "GeoFiles")

 #Loading the shapefiles
state_shape <- st_read("GeoFiles//cb_2018_us_state_500k.shp",
                       quiet = TRUE)

#calculating the rate of casualties 

population_full <-population_group|>
  group_by(state)|>
  summarize(total_population = sum(population))

total_casualties<-gun_data |>
  group_by(state)|>
  summarize(total_casualties = sum(n_injured+n_killed))


total_casualties<- total_casualties |>
  left_join(population_full)

total_casualties<-total_casualties |>
  mutate(casualty_rate = round((total_casualties/total_population)*100000))

merged_data <- left_join(state_shape, total_casualties,by = c("NAME"="state"))
state_shape <- st_transform(state_shape, "+proj=longlat +datum=WGS84")

# Define the bins and color palette
bins <- seq(0, 20, by = 3)
paletteNum <- colorBin(palette = "Reds", 
                       domain = merged_data$total_casualties, 
                       bins = bins)

labels <- merged_data |>
  dplyr::mutate(labels = paste0("<strong>", stringr::str_to_title(NAME), "</strong><br/>",
  "Casualty Rate: ", round(casualty_rate, 4)) |>
  lapply(htmltools::HTML)) |> 
  dplyr::pull(labels)

# Create Leaflet map
gun_violence_map <- leaflet()|>
  addProviderTiles(providers$Stadia.StamenTonerLite,
                   options = providerTileOptions(
    id = "mapbox.light")) |>
  setView(lng = -96.25, lat = 39.50, zoom = 3.5) |>
  addPolygons(
    data = state_shape,
    color = 'black',
    weight =  1,
    fillOpacity = 0.9,
    fillColor = ~paletteNum(merged_data$casualty_rate),
    label = labels,
    highlightOptions = highlightOptions(
      weight = 3,  
      color = "black",  
      fillOpacity = 1  
    )
  ) |>
  addLegend(
    pal = paletteNum,
    values = merged_data$total_casualties,
    title = "Total Casualties per 100k",
    position = "bottomleft"
  )

gun_violence_map